Bulletin of Mathematical Biology

Stratified epithelial tissues such as the skin epidermis maintain barrier integrity during development and homeostasis through the coordinated action of cell proliferation, differentiation, delamination, and tissue-scale mechanical forces. During development, the orientation of cell division within the basal layer plays a pivotal role in tissue stratification; however, the mechanical principles linking the orientation of the division plane to these processes across developmental stages remain poorly understood. Here, we expand a recently developed three-dimensional vertex model for stratified epithelia, composed of the basement membrane, basal, and suprabasal layers, to study the mechanical and structural impact of cell divisions with a wider range of orientations. The model integrates developmental stage via specific changes in heterotypic interfacial tensions (arising from actomyosin cortical contractility and adhesion molecules at the basal-suprabasal interface) and tissue stiffness that have been quantified previously in experiments. By systematically varying background mechanical parameters, we investigate how heterotypic tension, division orientation, and tissue fluidity collectively influence the outcome of cell division. Our goal is to uncover the strategies that the embryo may employ to generate stratified phenotypes at different developmental stages, recognizing that these strategies might evolve over time. Although our focus is on the embryonic developmental stages of the epidermis, this framework may also be extended to investigate transformed cells, such as in cancer, to explore how altered division orientation contributes to precancerous or transformed phenotypes.

Patterns in the vegetation across arid and semiarid regions may be explained as a form of self-organization driven by water scarcity, and are often modeled through reaction-diffusion dynamics. Recent work has shown that similar mathematical models generate patterns on networks. However, these studies have focused on idealized topologies with no reference to natural pattern-forming systems. Our study aims at bridging these two fields: we employ a physical reaction-diffusion vegetation model, and gradually modify the topology of the diffusion network by adding random shortcuts over a 2-dimensional grid, interpolating between a regular lattice and a random network. We found that network topology strongly shapes both the resulting vegetation patterns and the precipitation range that supports them. Three behavioral regimes emerge. On a regular lattice, high-regularity patterns develop reflecting local diffusion processes. On a random network, the system is dominated by global pressure towards homogenization yielding either a uniform state or a single patch. In the intermediate shortcut density range, as the network topology resembles a small world network, the interaction between the two scales of diffusion generates two kinds of disordered patterns: low-regularity patterns with a well-defined characteristic wavelength, and irregular patterns characterized by a broad patch size distribution. These disordered patterns resemble real-world observations and, in our model, they show different responses to changing precipitation. Although we focused on dryland vegetation, we suggest that network-mediated diffusion could lead to similar mechanisms in a wide variety of pattern-forming systems. HighlightsO_LIWe study vegetation pattern formation over different diffusion network topologies. C_LIO_LITwo kinds of stable disordered patterns states develop over small world topologies. C_LIO_LILow-regularity patterns with a well-defined characteristic wavelength. C_LIO_LIIrregular patterns characterized by a broad patch size distribution. C_LIO_LIThese different kinds of disordered states show different relations to precipitation. C_LI

Animals such as bees, ants, wasps, termites, and naked mole-rats live in colonies in which a single queen is the only female reproductive, an arrangement known as eusociality. Eusocial animals are known for their remarkably long lifespans. It has been argued that longevity becomes selected when queens are shielded from "external mortality". While such protection may contribute, we find a deeper reason: the eusocial reproduction strategy itself inherently creates selection for long lifespans. Lifespans typically reflect two processes: the baseline risk of death and the rate at which this risk increases with age. Each is a parameter in the Gompertz mortality equation. We show that the mathematical properties of eusocial reproduction lead to slowly-growing, older populations where selection acts more strongly on the rate at which risk increases than on the baseline risk. In addition, we show that channeling reproduction through a single female also selects for longevity, which we term the "queen effect". Thus, the dynamics of eusocial reproduction select for longer lifespan. More broadly, these results show that reproductive structure and population growth dynamics can fundamentally shape selection on lifespan, with implications outside eusocial systems as well.

The practical utility of many modern phylogenetic comparative methods can depend on how accurately mathematical models capture the evolutionary process of traits. Boucher and Demery (2016) described a new quantitative trait model, Brownian motion with reflective limits, that they anticipated might be of use in testing hypotheses about a particular sort of constraint on phenotypic character evolution. Since their analytic solution for the probability function under this bounded evolutionary scenario was not practical to evaluate for reasonably-sized trees, Boucher and Demery (2016) also identified a creative technique for computing the likelihood of their model. The basis of this methodology derives from the convergence of an equal-rates, symmetric, ordered Markov chain and continuous stochastic diffusion in the limit as the number of steps in our chain goes to {infty} (or, alternatively, as their widths decrease towards zero). We refer to this convergence in the limit as the discretized diffusion approximation or (more compactly) the discrete approximation. We realized that this discrete approximation of Boucher and Demery (2016) unlocked a number of additional models for the phylogenetic comparative analysis of discrete and continuous trait data, and we explore several of these in the present article. Specifically, we examine application of this discretized diffusion approximation to the threshold model from evolutionary quantitative genetics, to a new "semi-threshold" trait evolution model, to a joint model of discrete and continuous traits in which the discrete trait influences the rate of evolution of our continuous character, as well as a model where precisely the converse is true, and to a discrete character dependent multi-trend trended continuous trait evolution model. We conclude with some context for the origins of our article and discussion of other possible applications of this powerful approach.

During development, embryos store, transmit, and transform information to generate spatial patterns. Positional information (PI) quantifies how precisely cells form patterns at a given time, but cell motion has limited its application to static tissues. We introduce a framework for PI in dynamic tissues by decomposing mutual information between cells positions and properties over time into information flows contributing to PI preservation, loss and generation. These reveal information-theoretic signatures of ubiquitous developmental processes, including instruction, sorting and mixing, directly from data. Applying this framework to whole-embryo cell trajectories in Drosophila, mouse and zebrafish gastrulation, we provide local and global information-theoretic quantification of cell mixing and derive bounds on PI preservation imposed by tissue dynamics. Analyzing tissue flows as dynamical systems, we further show that morphogenesis structures mixing, preferentially preserving specific patterns. Finally, we derive inequality conditions for tracing generated PI to candidate information sources and distinguishing among alternative pattern-formation mechanisms, from programmed extracellular cues to self-organizing intercellular interactions.

We study an unexpectedly fast decay of motility in dense suspensions of Escherichia coli bacteria supplied with excess glucose under anaerobic conditions. The decrease in swimming speed occurs on a timescale inversely proportional to the cell concentration, and is associated with the secretion of organic acids by the bacteria. We show that the decay is driven by the progressive accumulation of non-ionised organic acids in the medium, and develop a chemical kinetic model that successfully predicts the swimming speed variations over a range of conditions in the presence of these acids. We further measure the internal pH of E. coli cells exposed to organic acids, and find that the speed decay coincides with sharp declines in internal pH and metabolic rate. Our findings identify an additional layer of motility control that can arise in complex environments even when motility genes are expressed and energy sources are abundant. This mechanism is likely relevant for understanding bacterial motility in habitats such as the human gut, where high densities of bacteria and organic acids are common.

PurposeColonic motility disorders, including diarrhea-predominant irritable bowel syndrome and slow-transit constipation, impose a major clinical burden. Although high-resolution colonic manometry reveals characteristic spatiotemporal motor patterns, such as high-amplitude propagating contractions and cyclic motor pattern in healthy individuals, these patterns are often altered or absent in disease. Understanding how these patterns arise from underlying pacemaker, neural, and mechanical mechanisms is essential for improving treatment strategies. MethodsWe developed a biophysical whole-colon model that integrates an Interstitial Cells of Cajal-inspired oscillator network, enteric nervous system reflexes, a pressure-gated modulation element motivated by rectosigmoid brake behavior, and a nonlinear tube law describing colon wall mechanics. The model simulates spatiotemporal pressure patterns along the colon and allows systematic variation of physiological parameters associated with pacemaker activity, neural reflex control, and distal gating. ResultsA small set of parameters reproduces three illustrative motility patterns corresponding to healthy motility, diarrhea-predominant irritable bowel syndrome, and slow-transit constipation. The simulated pressure maps recapitulate key features observed in high-resolution manometry, including propagation direction, regional patterning of contractions, and case-specific changes in amplitude and coordination. Sensitivity analysis suggests that proximal excitation strength and waveform morphology strongly influence global motility metrics. ConclusionOur study presents a simple, biophysical framework for reproducing clinically observed colonic motor patterns and exploring their disruption in disease. More broadly, the model may help interpret clinical manometry in mechanistic terms and support hypothesis-driven in silico studies of colonic motility disorders.

Functional spinal cord repair in zebrafish is governed by regeneration-favorable biochemical and mechanical cues within the lesion microenvironment. Alterations in extracellular matrix composition and stiffness are closely associated with axon regeneration. However, experimentally dissecting the interplay between mechanical signals and axonal regrowth in vivo remains technically challenging. Here, we present an agent-based modeling framework to simulate stiffness-mediated axonal growth trajectories across the lesion. We use this model to explore potential mechanisms underlying the characteristic growth patterns observed during zebrafish spinal cord regeneration. Computational predictions were qualitatively compared with confocal imaging data obtained from larval zebrafish. These phenomenological comparisons revealed a close agreement between simulated and experimentally observed axon growth, indicating that experimentally observed patterns could be governed by transient changes in the stiffness profile of the spinal cord and lesion microenvironment. Hence, our computational framework provides an in silico platform for investigating the role of mechanical cues in axon regeneration in the injured spinal cord.

A central goal of neuroscience is to understand how neural circuit architecture supports learning. While recent work has clarified the computational role of depth in sensory cortical hierarchies, it remains unclear why predominantly feedforward, non-convolutional circuits such as the cerebellum and olfactory system also contain multiple processing layers. Theoretical work in deep learning has shown that two-hidden-layer networks can achieve classification capacity that scales quadratically with the number of intermediate neurons, but these results rely on nonlocal synaptic optimization and are therefore difficult to reconcile with biological learning rules. Here, we show analytically and numerically that a two-hidden-layer network with feedforward gating can achieve quadratic capacity using local three-factor Hebbian learning when intermediate activity is sparse. This architecture supports efficient one-shot learning and, in settings where backpropagation requires many repeated weight updates, offers an advantage in learning speed. Beyond random perceptron tasks, the model also performs well on structured cerebellum-related tasks, including reinforcement-learning-based motor control. Mapping the model onto cerebellar microcircuitry further suggests functional roles for dendritic compartmentalization, branch-specific inhibition, and disinhibitory interneuron pathways. Together, these results extend the Marr-Albus-Ito framework by showing how the presence of multiple intermediate layers in cerebellum-like circuits can support fast, local, and high-capacity learning.

Abstract: Purpose: Recurrence in high grade glioma (HGG) predominantly occurs within the high dose radiation field, raising the question of whether treatment failure reflects limitations in radiation target delineation or is driven by intrinsic tumor biology. This study evaluated recurrence patterns following standard chemoradiotherapy and their treatment implications. Material and Methods: This retrospective single center study included 41 patients with histologically confirmed HGG treated with surgery followed by radiotherapy with concurrent and adjuvant temozolomide (TMZ). Patients were followed through August 2018; those with recurrence were included in the analysis. Recurrence patterns were classified based on their spatial relationship to the 60 Gy isodose line as central, infield, marginal, or distant. Survival outcomes were estimated using the Kaplan-Meier method and compared using the log rank test. Results: The most common pattern of recurrence was central (15 patients, 36.5%), followed by infield (11, 26.8%), distant (6, 14.6%), marginal (5, 12.1%), and multicentric (4, 9.8%). Central and in field recurrences (local failures) accounted for 26 patients (63%). Median overall survival (OS) was 27 months, and median progression-free survival (PFS) was 12 months. Survival differed significantly by recurrence pattern (log-rank p = 0.018), with marginal recurrence associated with more favorable outcomes. Conclusion: The predominance of central and infield recurrences within the high-dose region suggests that treatment failure in HGG is not solely explained by inadequate target delineation and may also be driven, in part, by intrinsic tumor biology, including radioresistant subpopulations and tumor heterogeneity. Future strategies may benefit from incorporating biologically guided approaches alongside optimization of radiation treatment parameters.

BackgroundLasers have wide applications in medicine and dermatology, but are associated with pain and anxiety, particularly in younger patients. Pain mitigation is often limited to topical anesthetics in the outpatient setting. Distraction techniques are limited by the need for ocular protection, which can include adhesive eye patches that can completely occlude vision. Virtual reality is effective at managing procedural pain and anxiety under other short medical procedures and is a promising tool for this population. ObjectiveThis trial aims to assess the safety, feasibility, and efficacy of Virtual Reality Pain Alleviation Therapeutic (VR-PAT) for pain management during outpatient laser procedures. Methods40 patients requiring outpatient laser therapy for at least two sessions will be recruited from a pediatric hospital in the midwestern United States for this crossover randomized, two-arm clinical trial with a 1:1 allocation ratio. During the first laser visit, the participant will be randomly assigned to either play the VR-PAT game during their procedure or wear the headset with a dark screen. Participants will answer questions about their pain (Numeric Rating Scale (NRS) 0-10), anxiety (State Trait Anxiety Inventory for Children, NRS 0-10, Modified Yale Preoperative Anxiety Scale (mYPAS)), and pain medication usage. Those playing the VR-PAT will additionally report simulator sickness symptoms and their experience playing the game. At their second laser visit, participants will crossover to the opposite intervention from their first visit. The primary outcomes are the difference in self-reported pain and anxiety between the two interventions. Feasibility outcomes include the proportion of screened patients who are eligible, consent, and complete both visits and adverse events reported. To evaluate the efficacy of pain reduction, composite scores of pain score, pain medication will be calculated for each laser visit. To evaluate the efficacy of anxiety reduction, the change of mYPAS scores will be compared between control and VR groups at each visit using Wilcoxon rank sum tests. All statistical analyses will follow the intention-to-treat principle in regard to intervention assignment at each visit. ResultsThe study was funded in January 2023 and began enrollment at that time. A total of n=44 participants were recruited and data collection was completed in November 2025, with n=40 subjects completing both visits. The sample was balanced with n=40 subjects using the intervention and participating in the control condition. The age range of the complete sample was 6 to 21 years at recruitment and was 55% female sex. Data analysis is in progress with final results planned for June 2026. ConclusionsFindings from this innovative randomized clinical trial will provide early evidence on the efficacy of the VR-PAT for reducing self-reported pain and anxiety during outpatient laser procedures. The results from this trial will inform a large-scale, multisite study. Trial RegistrationClinicalTrials.gov: NCT05645224 [https://clinicaltrials.gov/study/NCT05645224]

Target trial emulation (TTE) offers a principled way to estimate treatment effects using real-world observational data, but analyses of time-varying treatment strategies remain vulnerable to immortal time bias. The clone-censor-weight (CCW) approach is increasingly used to address this problem, yet key aspects of its causal interpretation and implementation remain unclear. In this work, we emulate a target trial using electronic health records (EHRs) to compare completion of a 3-dose 9-valent human papillomavirus vaccination (HPV) series within 12 months versus remaining partially vaccinated among vaccine initiators. We link CCW to the classic potential outcome framework in causal inference, evaluate the role of different weighting mechanisms, and account for within-subject correlation induced by cloning using cluster-robust variance estimation. Our study provides practical guidance for applying CCW in real-world comparative effectiveness studies to address immortal time bias and supports more rigorous and interpretable treatment effect estimation in TTE.

Mechanistic modeling has long been used as a tool to describe the dynamics of biological systems, especially cancer in response to treatment. Their key advantage lies in interpretability of relationships between input parameters and outcomes of interest. In contrast, machine learning techniques offer strong prediction performance, especially for high dimensional datasets that are common in oncology. Here, we employ a Mechanstic Learning framework that combines the advantages of both approaches by training machine learning models on mechanistic parameters inferred from clinical patient data. The mechanistic model (a Markov chain model) contains sixteen parameters that describe the rate of cell fate transitions that occur in patients with B-cell precursor acute lymphoblastic leukemia. The machine learning (a ridge logistic regression model) is trained on these parameters to predict two clinically-relevant features: BCR::ABL1 fusion gene status (positive or negative) and minimal residual disease status (positive or negative) post-induction chemotherapy. Model training is done in an iterative fashion to assess which (and how many) parameters are critical to maintain high predictive performance. Using machine learning models trained on the clinical flow-cytometry data, we find that the stem-like cell state alone is the most predictive feature for both BCR::ABL1-positive and MRD-positive disease, with combination scores (defined as the average of accuracy, balanced accuracy, and area under the curve) of 0.80 and 0.67, respectively. By comparison, mechanistic learning achieves comparable or improved combination scores for BCR::ABL1-positive and MRD-positive disease, with scores of 0.81 and 0.71, respectively, using only de-differentiation for BCR::ABL1 and primitive-state persistence together with differentiation-directed exit for MRD. Thus, the mechanistic-learning approach not only preserves predictive performance, but also provides a biological hypothesis for why stemness is predictive of these clinically relevant outcomes.

Audition is a fundamental sense, underlying critical human behaviors such as communication and recognition. Despite its importance, how the tuning and organization of receptive fields mature from childhood to adulthood in auditory cortex has not been directly measured. Through a gamified neuroimaging approach using functional MRI, we model population receptive field (pRF) tuning for frequency across human auditory cortex in both children and adults. In the same participants, we behaviorally quantify detection thresholds for different frequencies embedded in noise to understand how the functional development of human auditory cortex drives behavior. We find that while the tonotopic organization of pRFs is qualitatively present in early childhood, there is a protracted increase in the representation of low frequencies in tonotopic maps of primary auditory cortex. This maturation of pRF tuning appears to drive basic auditory behaviors, correlating with tone detection thresholds across participants. We also observe protracted development in secondary auditory regions, offering evidence for an anatomically-predictable tonotopic map posterior to Heschls Gyrus. These data provide a new avenue for studying the development of audition in the human brain and lay important groundwork for understanding atypical development in auditory processing disorders.

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has continuously evolved since its emergence in the human population in 2019. As of 1st August 2025, more than 1,700 Omicron subvariants have been designated by the Pango nomenclature system. The Pango nomenclature system designates a new lineage based on genetic and epidemiological information of SARS-CoV-2 strains. However, there is a possibility that strains that have similar genetic backgrounds and the same phenotype are given different Pango lineage names. In this paper, we propose a new algorithm, called FindPart-w, which can identify groups of viral lineages that share the same relative effective reproduction numbers. We introduced a new lineage replacement model, called the constrained RelRe model, which constrains groups of lineages to have the same relative effective reproduction numbers. The FindPart-w algorithm searches the equality constraints that minimise the Akaike Information Criterion of constrained RelRe models. Using hypothetical observation count data created by simulation, we found that the FindPart-w algorithm can identify groups of lineages having the same relative effective reproduction number in a practical computational time. Applying FindPart-w to actual real-world data of time-stamped lineage counts from the United States, we found that the Pango lineage nomenclature system may have given different lineage names to SARS-CoV-2 strains even if they have the same relative effective reproduction number and similar genetic backgrounds. In conclusion, this study showed that viruses that had the same relative effective reproduction number were identifiable from temporal count data of viral sequences. These findings will contribute to the future development of lineage designation systems that consider both genetic backgrounds and transmissibilities of lineages.

Invasive insect species pose serious threats to agriculture and ecosystems, with their spread increasingly accelerated by global trade and climate change. To support prevention and mitigation efforts, it is essential to map the regions where these pests can survive and thrive. Here, we apply MaxEnt, a leading species distribution modeling framework, to estimate current (2020) and future (2040-2060) suitable habitats for five major invasive insects across the contiguous United States: brown marmorated stink bug, corn earworm, spongy moth, root weevil, and spotted lanternfly. To account for an uncertain climatic future, these projections are generated under four shared socioeconomic pathways, which reflect a range of plausible climate change scenarios. Beyond forecasting distributions, we examine several key modeling decisions, especially those often overlooked in practice. In particular, we find that background sampling strategies play a critical role in model calibration and that a hybrid sampling approach with a moderate buffer bias provides better predictive accuracy. We also show that permutation importance scores, commonly used to rank environmental variables, are highly sensitive to small changes in the background data and should be interpreted with caution. Finally, to bridge the gap between ecological modeling and applied machine learning, we provide a self-contained, math-focused background to MaxEnt aimed at practitioners outside of traditional ecological fields. Overall, this work delivers reproducible modeling workflows and critical insights into building robust, transparent, and ecologically meaningful MaxEnt models for climate-informed species distribution analysis.

Lesbian, gay, bisexual, transgender, queer, intersex, and asexual (LGBTQIA+) students continue to face violence, exclusion, and barriers at school, including in STEM education. A key underexamined factor in diversity, equity, and inclusion (DEI) efforts is the content of the life science curriculum, which is uniquely positioned to reinforce or refute bioessentialist, binary, and heteronormative biases. Outdated science curricula not only conflict with current scientific evidence but can also perpetuate beliefs that contribute to sexism and LGBTQIA+ marginalization. To address this, we designed four gender and sexual diversity (GSD)-inclusive biology activities, aligned with NGSS standards, and informed by inclusive curriculum frameworks. Using a mixed-methods approach, we studied 127 high school students who participated in two or more inclusive biology activities. Surveys conducted before and after implementation showed significant reductions in essential, binary beliefs about sex and gender, and increases in affirming attitudes toward sex and gender diversity. Interviews conducted after implementation further revealed differences between LGBTQIA+ and straight students conceptualizations of biological sex. Our findings demonstrate that even brief curriculum interventions can shift student attitudes, although we hope future studies will explore the impact of sustained interventions. Updating life science instruction is essential for educational equity and scientific accuracy.

Background and AimsFenofibrate is widely prescribed for hyperlipidaemia and has been associated with rare but severe cases of drug-induced liver injury (DILI), yet its effects on liver sinusoidal endothelial cells (LSECs) remain to be investigated. LSECs maintain a highly permeable specialized sinusoidal barrier characterized by transcellular pores (fenestrations), regulating the bidirectional transfer of circulating compounds to and from the hepatocytes. As drug-induced alterations in fenestration architecture could influence xenobiotic access to hepatocytes, these changes may modulate pathways associated with DILI. Understanding the effects of fenofibrate on LSEC ultrastructure may therefore provide insights into previously underexplored endothelial contributions to hepatic drug responses. MethodsBoth fenofibrate and its active metabolite, fenofibric acid, were evaluated for their effects on LSEC ultrastructure, mechanical properties, and functional markers. Atomic force microscopy (AFM) and scanning electron microscopy (SEM) and were used to quantify fenestration architecture. AFM was additionally used to measure cellular mechanical properties, which were interpreted in the context of fluorescence-based quantification of cytoskeletal organization. Gene expression, viability, and cytotoxicity were assessed using PCR-based and biochemical assays. ResultsFenofibrate reduced fenestration number and porosity at both tested concentration (10, and 25 {micro}M). It also decreased the apparent Youngs modulus of LSECs, accompanied by changes in tubulin and actin architecture, without detectable cytotoxicity. In contrast, treatment with fenofibric acid did not result in significant structural or mechanical effects on LSECs, even at higher concentrations. ConclusionsTogether, these data identify LSECs as a drug-responsive hepatic cell type for fenofibrate, suggesting that LSECs could represent an underrecognized contributor to the complex, multifactorial processes underlying DILI. This work provides a framework for evaluating endothelial contributions to fenofibrate-associated liver effects in more complex models. O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=105 SRC="FIGDIR/small/718907v1_ufig1.gif" ALT="Figure 1"> View larger version (51K): org.highwire.dtl.DTLVardef@1d3f60corg.highwire.dtl.DTLVardef@bea13aorg.highwire.dtl.DTLVardef@14b27d8org.highwire.dtl.DTLVardef@124e0d3_HPS_FORMAT_FIGEXP M_FIG Fenofibrate reduces LSEC fenestrations and metabolic activity at higher concentrations, while its metabolite, fenofibric acid, does not affect LSEC, regardless of its concentration. C_FIG

Widespread genomic sequencing efforts have characterized the molecular foundations of the different cancers. By combining these genomic data in a manner proportional to the population-level abundances of these different cancers, we estimate the overall abundances of each observed missense and nonsense mutation within the U.S. cancer patient population. We find BRAF V600E (5.2%) is the most common mutation in the cancer patient population, TP53 R175H (1.5%) is the most common tumor suppressor mutation, and APC R876X (0.4%) is the most common nonsense mutation. These values differ largely and significantly from what would be found in a typical pan-cancer analysis, where different cancer types are included out of proportion to population level incidence. We present the full ordered lists of population-level abundances for specific missense and nonsense mutations, and we demonstrate the value of these data by further analyzing high priority genes (e.g., TP53, KRAS, BRAF) and pathways (e.g., RTK/RAS, PI3K, and WNT/{beta}-catenin). Overall, this information is a resource that should benefit the basic science, translational, and clinical cancer research communities.

Pathogens that persist subclinically across many wildlife populations can drive mass mortality in others. Mass mortality is often abrupt, and the timing can be difficult to predict from host or habitat features alone. In a recent field study tracking ranavirus epizootics in wood frog (Rana sylvatica) breeding ponds, we found that no environmental or biotic feature reliably predicted die-off occurrence or timing. Instead, the trajectory of viral accumulation in the water column was the strongest dynamic predictor of mass mortality. Infected hosts shed virus throughout epizootics, but the influence of waterborne viral concentration on disease progression was apparent only near die-off onset. This pattern suggests a potential threshold-dependent feedback operating through the shared viral environment. Here, we develop a compartmental model linking waterborne viral concentration to the rate at which subclinical infections progress to clinical, high-shedding states within already-infected hosts. We show that a dose-dependent progression model generates the two-phase epizootic trajectory observed in natural die-offs: prolonged subclinical circulation followed by abrupt clinical transition after environmental virus crosses an escalation threshold. The model exhibits a sharp phase transition between subclinical circulation and mass mortality, governed mainly by the clinical-to-subclinical shedding ratio, host density, and pond volume. Existing explanations for die-off variation emphasize individual-level susceptibility, but our model demonstrates that dose-dependent environmental feedback, a mechanism not previously formalized at the population level, can generate the transition from subclinical infection to mass mortality without invoking individual variation in host susceptibility. This mechanism may apply in any system where hosts share a bounded environment, pathogen dose influences disease severity, and pathogen shedding increases with disease progression.